function output = momentsModelEval(yieldData,shortRate,shortRateOneYear,econAct,...
    ygap,thresholdLevel,maturities,sampletime_Per_Maturitytime,bootOn)

regime_1           = real(econAct) >= thresholdLevel; %a boom
regime_2           = real(econAct) < thresholdLevel;  %a recression
meanYield          = mean(yieldData,1);
meanYield_regime1  = mean(yieldData.*repmat(regime_1,1,size(yieldData,2)),1);
meanYield_regime2  = mean(yieldData.*repmat(regime_2,1,size(yieldData,2)),1);
stdYield           = std(yieldData,0,1);
stdYield_regime1   = std(yieldData.*repmat(regime_1,1,size(yieldData,2)),0,1);
stdYield_regime2   = std(yieldData.*repmat(regime_2,1,size(yieldData,2)),0,1);

if bootOn == 1
    numBoot = 1000;
    numYields = size(yieldData,2);
    meanYield_bs         = zeros(numYields,numBoot);
    meanYield_regime1_bs = zeros(numYields,numBoot);
    meanYield_regime2_bs = zeros(numYields,numBoot);
    stdYield_bs          = zeros(numYields,numBoot);
    stdYield_regime1_bs  = zeros(numYields,numBoot);
    stdYield_regime2_bs  = zeros(numYields,numBoot);
    bsdata=block_bootstrap([econAct yieldData],12*10,numBoot);
    for s=1:numBoot
         econ_bs      = bsdata(:,1,s);
         yields_bs    = bsdata(:,2:end,s);
         regime_1_bs  = real(econ_bs) >= thresholdLevel; %a boom
         regime_2_bs  = real(econ_bs) < thresholdLevel;  %a recression
         meanYield_bs(:,s)          = mean(yields_bs,1);
         meanYield_regime1_bs(:,s)  = mean(yields_bs.*repmat(regime_1_bs,1,size(yields_bs,2)),1);
         meanYield_regime2_bs(:,s)  = mean(yields_bs.*repmat(regime_2_bs,1,size(yields_bs,2)),1);
         stdYield_bs(:,s)           = std(yields_bs,0,1);
         stdYield_regime1_bs(:,s)   = std(yields_bs.*repmat(regime_1_bs,1,size(yields_bs,2)),0,1);
         stdYield_regime2_bs(:,s)   = std(yields_bs.*repmat(regime_2_bs,1,size(yields_bs,2)),0,1);
    end
    meanYield_bootSE         = std(meanYield_bs,0,2);
    meanYield_regime1_bootSE = std(meanYield_regime1_bs,0,2);
    meanYield_regime2_bootSE = std(meanYield_regime2_bs,0,2);
    stdYield_bootSE          = std(stdYield_bs,0,2);
    stdYield_regime1_bootSE  = std(stdYield_regime1_bs,0,2);
    stdYield_regime2_bootSE  = std(stdYield_regime2_bs,0,2);
end

LPYI               = CampbellSchillerRegression(yieldData,shortRate,zeros(size(yieldData)),zeros(size(yieldData)),maturities,sampletime_Per_Maturitytime);
LPYI_Regime        = CampbellSchillerRegression_Regime(yieldData,shortRate,econAct,thresholdLevel,maturities,sampletime_Per_Maturitytime);
[xhrOnSlope,xhrTmp,termSpreads]= forecastRegressionYieldSpread(yieldData,shortRate,shortRateOneYear,maturities,regime_2,sampletime_Per_Maturitytime);
AR1_regime         = AR1_regimeSwitching(yieldData,maturities,econAct,thresholdLevel);

% Regressing realized returns on changes in the short rate
returnsOnShortRate = regressReturnsShortRate(yieldData,shortRateOneYear,maturities,regime_2,sampletime_Per_Maturitytime);

% Forecasting econAct using slope
pmiOnSlope     = regressMacroOnSlope(yieldData,shortRateOneYear,maturities,regime_2,econAct,sampletime_Per_Maturitytime);
pmiDummyOnSlope= regressMacroOnSlope(yieldData,shortRateOneYear,maturities,regime_2,regime_2,sampletime_Per_Maturitytime);
ygapOnSlope    = regressYGapOnSlope(yieldData,shortRateOneYear,maturities,ygap,sampletime_Per_Maturitytime);

% Finding the dates where the recession starts
LastObsBeforeRec = [(regime_2(1:end-1) == 0 & regime_2(2:end) == 1)*1;0];
% To address the case where we start in a recession
if regime_2(1) == 1 && regime_2(2) == 1
    LastObsBeforeRec(1) = 1;
end
endRec           = [(regime_2(1:end-1) == 1 & regime_2(2:end) == 0)*1;0];
endRec(1,1)      = 0;
% To address the case where we end the sample in a recession
if regime_2(end) == 1 
    endRec(end) = 1;
end
dataDyieldRec    = yieldData(endRec==1,:)-yieldData(LastObsBeforeRec==1,:);
avgdYieldsRec    = mean(yieldData(endRec==1,:)-yieldData(LastObsBeforeRec==1,:),1);
idx              = 1:length(regime_2);
recDuration      = idx(endRec==1)-idx(LastObsBeforeRec==1);

xhr           = [nan(sampletime_Per_Maturitytime,size(xhrTmp,2));xhrTmp];
[T,maxMat] = size(xhr);
numRecPer = 8;
xhrAfterRecStart        = zeros(sum(LastObsBeforeRec),maxMat,numRecPer);
xhrFitAfterRecStart     = zeros(sum(LastObsBeforeRec),maxMat,numRecPer);
termSpreadAfterRecStart = zeros(sum(LastObsBeforeRec),maxMat,numRecPer);
DyieldsAfterRecStart    = zeros(sum(LastObsBeforeRec),maxMat,numRecPer);
idx = 0;
for t=1:T
    if LastObsBeforeRec(t,1) == 1
        idx = idx + 1;
        for i=1:numRecPer
            if t+i <= T
                xhrAfterRecStart(idx,:,i)        = xhr(t+i,:);
                xhrFitAfterRecStart(idx,:,i)     = xhrOnSlope.betaRegime(2,:)+termSpreads(t+i,:).*xhrOnSlope.betaRegime(4,:);
                termSpreadAfterRecStart(idx,:,i) = termSpreads(t+i,:);
                DyieldsAfterRecStart(idx,:,i)     = yieldData(t+i,:)-yieldData(t,:);
            else
                xhrAfterRecStart(idx,:,i)        = xhr(T,:);
                xhrFitAfterRecStart(idx,:,i)     = xhrOnSlope.betaRegime(2,:)+termSpreads(T,:).*xhrOnSlope.betaRegime(4,:);
                termSpreadAfterRecStart(idx,:,i) = termSpreads(T,:);
                DyieldsAfterRecStart(idx,:,i)     = yieldData(T,:)-yieldData(t,:);
            end
        end
    end
end

meanXhrAfterRecStart        = zeros(maxMat,numRecPer);
meanXhrAfterRecStartSE      = zeros(maxMat,numRecPer);
meanXhrFitAfterRecStart     = zeros(maxMat,numRecPer);
meanXhrFitAfterRecStartSE   = zeros(maxMat,numRecPer);
meanTermSpreadAfterRecStart = zeros(maxMat,numRecPer);
meanTermSpreadAfterRecStartSE = zeros(maxMat,numRecPer);
meanDyieldsAfterRecStart    = zeros(maxMat,numRecPer);
meanDyieldsAfterRecStartSE  = zeros(maxMat,numRecPer);

for k1=1:maxMat
    for k2=1:numRecPer
        res = nwest(xhrAfterRecStart(:,k1,k2),ones(sum(LastObsBeforeRec),1),numRecPer);
        meanXhrAfterRecStart(k1,k2) = res.beta;
        meanXhrAfterRecStartSE(k1,k2) = res.beta/res.tstat;
        
        res = nwest(xhrFitAfterRecStart(:,k1,k2),ones(sum(LastObsBeforeRec),1),numRecPer);
        meanXhrFitAfterRecStart(k1,k2) = res.beta;
        meanXhrFitAfterRecStartSE(k1,k2) = res.beta/res.tstat;
        
        res = nwest(termSpreadAfterRecStart(:,k1,k2),ones(sum(LastObsBeforeRec),1),numRecPer);
        meanTermSpreadAfterRecStart(k1,k2) = res.beta;
        meanTermSpreadAfterRecStartSE(k1,k2) = res.beta/res.tstat;
        
        res = nwest(DyieldsAfterRecStart(:,k1,k2),ones(sum(LastObsBeforeRec),1),numRecPer);
        meanDyieldsAfterRecStart(k1,k2) = res.beta;
        meanDyieldsAfterRecStartSE(k1,k2) = res.beta/res.tstat;
    end
end




%% saving output
output.meanYield         = meanYield;
output.meanYield_regime1 = meanYield_regime1;
output.meanYield_regime2 = meanYield_regime2;
output.stdYield          = stdYield;
output.stdYield_regime1  = stdYield_regime1;
output.stdYield_regime2  = stdYield_regime2;
if bootOn == 1
    output.meanYield_bootSE         = meanYield_bootSE';
    output.meanYield_regime1_bootSE = meanYield_regime1_bootSE';
    output.meanYield_regime2_bootSE = meanYield_regime2_bootSE';
    output.stdYield_bootSE          = stdYield_bootSE';
    output.stdYield_regime1_bootSE  = stdYield_regime1_bootSE';
    output.stdYield_regime2_bootSE  = stdYield_regime2_bootSE';
end
output.LPYI              = LPYI;
output.LPYI_Regime       = LPYI_Regime;
output.AR1_regime        = AR1_regime;
output.xhrOnSlope        = xhrOnSlope;
output.returnsOnShortRate= returnsOnShortRate;
output.pmiOnSlope        = pmiOnSlope;
output.pmiDummyOnSlope   = pmiDummyOnSlope;
output.ygapOnSlope       = ygapOnSlope;
output.maturities        = maturities;
output.avgdYieldsRec     = avgdYieldsRec;
output.dataDyieldRec     = dataDyieldRec;

output.meanXhrAfterRecStart        = meanXhrAfterRecStart;
output.meanXhrAfterRecStartSE      = meanXhrAfterRecStartSE;
output.meanXhrFitAfterRecStart     = meanXhrFitAfterRecStart;
output.meanXhrFitAfterRecStartSE   = meanXhrFitAfterRecStartSE;
output.meanTermSpreadAfterRecStart = meanTermSpreadAfterRecStart;
output.meanTermSpreadAfterRecStartSE = meanTermSpreadAfterRecStartSE;
output.meanDyieldsAfterRecStart    = meanDyieldsAfterRecStart;
output.meanDyieldsAfterRecStartSE  = meanDyieldsAfterRecStartSE;
output.recDuration                 = recDuration;

end